Spatial coupling of particle and fluid models for streamers: where nonlocality matters 
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Particle models for streamer ionization fronts contain correct electron energy distributions, run- 
away effects and single electron statistics. Conventional fluid models are computationally much more 
efficient for large particle numbers, but create too low ionization densities in high fields. To combine 
their respective advantages, we here show how to couple both models in space. We confirm that 
the discrepancies between particle and fluid fronts arise from the steep electron density gradients in 
the leading edge of the fronts. We find the optimal position for the interface between models that 
minimizes computational effort and reproduces the results of a pure particle model. 

PACS numbers: 52.80.-s, 52.80.Mg, 52.65.Kj 



Streamers generically occur in the initial electric break- 
down of long gaps. They are growing filaments of weakly 
ionized nonstationary plasma; they are produced by a 
sharp ionization front that propagates into non-ionized 
matter within a self-enhanced electric field. Streamers 
are used in industrial applications such as lighting 1], 
gas and water purification or combustion control 4|, 
and theyoccur in natural processes as well such as light- 
ning [5J, [6| or transient luminous events in the upper at- 
mosphere @, Important recent questions concern (i) 
propagation and branching of streamers @ and the role 
of avalanches created by single electrons, (n) the elec- 
tron energy distribution in the streamer head and the 
subsequent gas chemistry that is used in the above appli- 
cations, as well as (in) runaway electrons and X-ray gen- 
eration, possibly in the streamer zone of lightning lead- 
ers [1, fiof . The present paper deals with the efficient 
simulation of these problems. 

Monte Carlo particle simulations [Il|, [Tl] model these 
effects as they contain the full microscopic physics; the 
deterministic electron motion between collisions is calcu- 
lated and collisions of electrons with neutrals are treated 
through a Monte Carlo process with appropriate statis- 
tical weights. The particle model includes the complete 
electron velocity and energy distribution as well as the 
discrete nature of particles. However, a drawback of such 
models is that the required computation resources grow 
with the number of particles and eventually exceed the 
CPU space of any computer. This difficulty is counter- 
acted by using superparticles carrying the charge and the 
mass of many physical particles, but superparticles in 
turn create unphysical fluctuations and stochastic heat- 
ing 0. 

Streamers are therefore mostly modeled as fluids (see 
e -g- [HL EE EIj HzL Ell) smce a fluid model is computa- 
tionally much more efficient. In the case of a negative 
discharge in a pure non-attaching gas like nitrogen, it 
consists of continuity equations for the densities of elec- 
trons n e and positive ions n p coupled to the Poisson equa- 
tion for the electric field E. The electron mobility \x 
and diffusion matrix D and the impact ionization rate a 
are calculated from microscopic scattering and transport 
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Figure f : (Color online) The streamer ionization front, that 
here is indicated by the electron density n e (z), and its presen- 
tation by particle or fluid model in different spatial regions. 



models like the Boltzmann equation [19| or directly from 
Monte Carlo simulations as, e.g., in [2Cj. In streamer 
calculations, it is generally assumed that these transport 
and reaction coefficients are functions of the local electric 
field. 

We recently have compared the properties of streamer 
ionization fronts of particle models and conventional fluid 
models [20j for negative planar fronts in nitrogen; the 
transport coefficients for the fluid model were generated 
from swarm experiments in the particle model. We found 
that the models agree reasonably for fields up to 50 
kV/cm at standard temperature and pressure, but that 
differences increase with increasing electric field. For ex- 
ample, in a field of 200 kV/cm, the ionization level behind 
the front is 60% higher in the particle model than in the 
fluid model. We have related this to the fact that the 
electron energies and, consecutively, the ionization rates 
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in the leading edge of the front are considerably higher 
in the particle than in the fluid model; they are actu- 
ally at the edge of runaway. We found that this effect 
is due to the strong density gradients in the front, and 
not due to field gradients. So for high fields and con- 
secutively strong density gradients at the streamer tip, 
there is a clear need for particle simulations, and par- 
ticles, rather than superparticles, should be used to get 
physically realistic density fluctuations when modeling, 
e.g., the branching process of a streamer. 

The basic idea of the present paper is demonstrated 
in Fig. [TJ namely to follow the single electron dynamics 
in the high field region of the streamer where the elec- 
tron density gradient is steep, and to present the interior 
region with large numbers of slower electrons through a 
fluid model with appropriate transport coefficients. As 
in [2(|, we study negative streamers in nitrogen, and we 
simplify the notation by refering to standard temperature 
and pressure though the model trivially scales with gas 
density. The particle and the fluid model by themselves 
are taken as described in detail in [2(|. But how should 
particle and fluid model be coupled in space? And where 
should the interface between the models be located to get 
fast, but reliable results? The answers to these questions 
will be given below. They required us to identify cor- 
rectly the spatial region where particle and fluid model 
deviate, and this allowed us to then compute the full 
electron physics efficiently in the relevant region. 

When coupling the models, the model interface should 
move with the ionization front; this keeps the total num- 
ber of electrons limited and superparticles need not be 
introduced. The position of the interface can be cho- 
sen either according to the electron densities or to the 
electric field. As the electron densities fluctuate stronger 
than the electric field, we relate the position of the model 
interface to the electric field. More precisely, the inter- 
face is placed where the local field E is a given fraction x 
of the maximal field E + : -©interface = x E + . By varying 
x, the region modeled by particles can be varied. 

To properly handle the interaction of two models, we 
introduced a so called "buffer region" where a parti- 
cle model coexists with a fluid model. The separation 
of the full computational region into fluid, particle and 
buffer region is indicated in Fig , f . Buffer regions have 
been introduced in [2l|, [2^, UK for rarefied gases cou- 
pling a direct simulation Monte Carlo (DSMC) scheme 
to the Navier-Stokes equations, and in other applica- 
tions [24], Physical observables are evaluated from 
the fluid model in its whole definition region up to the 
model interface. Beyond that point, the particle model 
is used. The particle model extends back beyond the 
model interface into the buffer region where particle and 
fluid model coexist; it supplies particle fluxes to the fluid 
model on the model interface. However, correct particle 
fluxes require correct particle statistics within the buffer 
region whose length should be as small as possible to re- 
duce computation costs, but larger than the electron en- 
ergy relaxation length [20(. In many cases, new particles 




Figure 2: Flow chart for one time step from t n to t n +i in the 
complete hybrid calculation. 



need to be introduced into the buffer region, that have 
to be drawn from appropriate distributions in configura- 
tion space. This would pose a particular problem since a 
Boltzmann or even a Druyvesteyn distribution can be in- 
accurate. But for negative streamers, where electrons on 
average move somewhat slower than the ionization front, 
the electron loss at the end of a sufficiently long buffer 
region does not affect the calculation of particle fluxes 
at the model interface. Therefore the particle loss at the 
end of the buffer region can be ignored and new electrons 
do not need be created artificially. 

In more detail, the calculation is performed as follows. 
One hybrid computation step from t n to t n +i is described 
in the flow chart in Fig. [2] The electric field E, the elec- 
tron and ion densities n e and n p in the fluid region and 
the kinetic information of particles in the particle and 
buffer region are given at time step t n . First, the po- 
sitions and velocities of all old and the newly generated 
particles are updated to time step i n +i in the particle and 
in the buffer region. Their collisions during this time step 
are treated stochastically and their new velocities and po- 
sitions are calculated by solving the equation of motion. 
The number of electrons crossing the model interface dur- 
ing this time step is recorded. This particle flux across 
the interface provides the required boundary condition 
for calculating the evolution of the densities in the fluid 
region up to the same time t n+ \. The particles are at 
arbitrary positions, but densities and field are calculated 
here on the same numerical grid. Therefore the particles 
in the particle region are averaged to densities on the nu- 
merical grid, and then the electric field at time is 
calculated from the Poisson equation everywhere. (The 
charge density in the buffer region is taken from the fluid 
model, and the particles in the buffer region only serve 
to generate correct fluxes for the fluid region.) Finally, 
the position of the model interface is updated to its new 
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Figure 3: Electron density n e and ion density n v (solid and 
dotted line) and electric field strength \E\ (dash-dotted line) 
in a streamer ionization front in nitrogen in a field of E + = 
— 100 kV/cm at standard temperature and pressure within a 
pure particle model. Our units are related to other commonly 
used units like 1 kV cm" 1 bar" 1 = 1.316 V cm" 1 torr" 1 = 
0.424 Townsend at T=300 K. Below we will model the leading 
edge of the front by a particle model and the streamer interior 
by a fluid model where the model interfaces are located at 
^interface = x E+ with x = 0.6, 0.7, 0.8, 0.9, and 0.98. These 
interface positions for Figs. 2 and 3 are marked by vertical 
lines, with "o" for the fields and with "+" for the densities. 



position at time fn+ij it can stay where it was or move 
one grid size forward. The buffer region moves with it. 
All particles that now are neither in the particle nor in 
the buffer region, are removed from the particle list. 

In the particle model, a standard PIC/MCC (Particle 
in Cell/Monte Carlo collision) method is implemented. 
At each time step of length At = 0.3 ps, particles in the 
particle regions are mapped to densities on a uniform 
grid with mesh size Al = 2.3 /im. Meanwhile, the fluid 
equations are solved in the fluid region of the same grid; 
discretization and grid dep endence of the fluid model are 
discussed in detail in [14j |. The charge density n p — n e 
then can be obtained everywhere and the electric field is 
calculated on this grid. The size of the time step and the 
grid size are chosen such that the ionization front need 
several time steps to move over one A£, e.g., 10 At at 
100 kV/cm and 3 At at 200 kV/cm. 

The length of the buffer region is another crucial factor 
in the hybrid computation. A buffer region with length 
of 32 A£ has been used in the present simulations, which 
ensure a reliable flux around the model interface and sta- 
ble results of hybrid simulations. The length of the buffer 
region is much larger than the energy relaxation length 
found in [2(j. The long buffer region does not bring a 
heavy burden to the simulation of the planar front sys- 
tem but will considerably reduce the computational effi- 
ciency in a more complex geometry. Therefore, the min- 
imal length of the buffer region as well as other features 
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Figure 4: Electron and ion densities in the coupled model 
(thick lines) in a field of E + = -100 kV/cm with model 
interfaces at interface = x E + with (a) x — 0.6, (b) x = 
0.9, and (c) x — 0.98; these interface positions are indicated 
by vertical dashed lines. The densities in the particle model 
(electrons: solid, ions: dot-dashed) and in the fluid model 
(electrons: dashed, ions: dotted) are shown as well; they are 
discussed in 12011. 



of fluid and particle models shall be investigated in more 
detail in a furture paper. 

We first show the simulation results of this coupled 
model for a front propagating into a field of E + = —100 
kV/cm, and with the model interface located at x = 0.6, 
0.9, and 0.98; the positions of these interfaces are indi- 
cated in Fig. [3l The field ahead of the front is fixed, and 
the system is always taken long enough that effects at 
the outer boundaries are not felt. The coupled model 
generates different electron and ion densities behind the 
ionization front as shown in Fig. SJ for x = 0.6, the den- 
sity is as in the particle model; for x = 0.98, it is as in the 
fluid model; and for x = 0.9, it takes some intermediate 
value. We conclude that the solution of the pure particle 
model can be replaced by the coupled model, if a suffi- 
ciently large region of the ionization front with its steep 
gradients is covered by the particle model, and that the 
coupling to the fluid model behind that region does not 
cause numerical artifacts. This confirms the discussion 
in [20j : it is indeed the high electron density gradient 
that causes an electron energy overshoot and a higher 
ionization rate in the leading edge of the particle front. 
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Figure 5: Relative density difference behind the front 
( n e, P art ~ n e, coup) / n e, part between the particle model and the 
coupled model as a function of the applied electric field E + 
and of the position of the model interface -Einterfacc = x ■ E + 
with x = 0.6 (□), 0.7 (x), 0.8 (■), 0.9 (*), 0.98 (o) and 1.0 
(+). The last case corresponds to the fluid model. (Note that 
a density difference of 60% relative to the fluid model corre- 
sponds to a density difference of 37.5% relative to the particle 
model.) 



The coupled model also confirms that the field gradients 
do not play a role in causing the density discrepancy be- 
tween the fluid and the particle model as the field keeps 



varying across the model interface in the coupled model. 

Having analyzed the ionization front propagating into 
a field of E + = — 100 kV/cm, we now summarize the 
results for fields ranging from —50 to —200 kV/cm. 
Fig. [5] shows the discrepancy between particle and cou- 
pled model on the most sensitive observable [20l ] , namely 
on the relative difference (n~ part - n~ coup ) /n- part of 
the saturated electron density n~ behind the ionization 
front. This quantity is shown as a function of the electric 
field E + and of the position of the model interface pa- 
rameterized again by x. The figure shows that for higher 
E + , the parameter x needs to be smaller. This shift of re- 
quired interface position relative to E + corresponds to a 
shift of the maximal electron density relative to E: both 
for E+ = -50 kV/cm and for E+ = -200 kV/cm, the 
particle and the coupled model agree well, if the model 
interface lies at the maximum of the electron density and 
therefore covers the complete steep gradient region; this 
is the case at E = 0.8 E+ for E+ = -50 kV/cm and at 
E = 0.35 E+ for E+ = -200 kV/cm. 

Coupling particle and fluid models in space with 
varying interface positions confirms our prediction (20| 
that the density discrepancies between particle and fluid 
model are due to the strong density gradients in the lead- 
ing edge of the front. This investigation also lays the 
basis for constructing a fully 3D coupled particle-fluid 
model where the fields ahead the ionization front are 
changing in space and time. 
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